Electron-phonon coupling as an order-one problem 
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The coupling between electrons and phonons plays important roles in physics, chemistry and 
biology. However, the accurate calculation of the electron-phonon coupling constants is computa- 
tionally expensive as it involves solving the Schrodinger equation for 0{N) nuclear configurations, 
where is the number of nuclei. Herein we show that by considering the forces on the nuclei 
caused by the addition or subtraction of an arbitrarily small electronic charge one may calculate 
the electron-phonon coupling constants from 0(1) solutions of the Schrodinger equation. We show 
that Janak's theorem means that this procedure is exact within the density functional formalism. 
We demonstrate that the 0(1) approach produces numerically accurate results by calculating the 
electron-phonon coupling constants for a series of molecules ranging in size from H2 to Ceo- We 
use our approach to introduce a computationally fast approximation for the adiabatic ionisation 
potentials and electron affinities which is shown to be accurate for large molecules. We also show 
that our approach allows for the calculation isotope effects in 0(0) time and discuss the deuteration 
driven Mott transition in K-(BEDT-TTF)2Cu[N(CN)2]Br. 

PACS numbers: 



Electrons in condensed matter feel two forces: a di- 
rect Coulomb interaction between themselves and an 
Columbic interaction with the atomic nuclei. As the 
vibrations of the nuclei are quantised the interaction 
with the lattice takes the form of the electron-phonon 
(or electron-vibron) interaction. Electron-phonon inter- 
actions plays important roles across science from physics 
(electron-phonon interactions can give rise to supercon- 
ductivity 1', spin and charge density waves 2], polaron 
formation 3] and piezoelectricity 0]), through chemistry 
(such as in electron transfer processes , Jahn- Teller ef- 
fects 0, spectroscopy , stereochemistry activation 
of chemical reactions j^l and catalysis 5]) to biology (for 
example electron-phonon interactions play an important 
role in photoprotection 1^, photosynthesis fi\ and vision 
0). It is therefore clear that one of the central tasks for 
condensed matter theory and theoretical chemistry is to 
accurately calculate electron-phonon coupling constants. 
However, until now, this has been a computationally ex- 
pensive task. 

The expense of calculating the electron-phonon cou- 
pling constants arises from the large number of phononic 
modes that need to be considered. For example, a (non- 
linear) molecule with atoms has M — 3N — 6 phononic 
modes. Within the Born-Oppenheimer approximation, 
the standard, frozen phonon, approach to calculating 
the electron-phonon coupling constants requires the elec- 
tronic eigenstates to be calculated for M nuclear config- 



urations. Below we demonstrate how the number of cal- 
culations can be reduced from M to 1 without increasing 
the difficulty of the electronic problem. 

Let us begin by reviewing the basic problem of 
electron-phonon coupling in the semiclassical limit. In 
this limit the phonons are represented by harmonic oscil- 
lators, thus the Hamiltonian is 
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where pi is the momentum of the z*'* mode, nii is the 
effective mass of the i*^ mode, ki is the spring constant 
of the z"* harmonic oscillator, Xi is the displacement of 
the z"* mode from its equilibrium position, the gi are 
the electron-phonon coupling constants, n is the number 
of electrons and 7^° is the Hamiltonian of the electrons 
in the absence of phonons. In what follows it will be 
useful to introduce dimensionless normal coordinates for 
the phonons given by Qi = Xi^fm^pJifh,, where LOi = 
y^ki/mi is the frequency of the phonon, dimensionless 
momenta. Pi — pi y^rrii/huji, and dimensionless electron- 
phonon coupling constants, gi ~ gi/ y^2mihujf . Thus, 
the Hamiltonian may be written as 



(2) 



We may quantise this Hamiltonian by introducing the 



Bosonic operators a^p which annihilate (create) a phonon 
on the in the z*'' mode and the Fermionic operators 
which annihilate (create) an electron in the state ^ such 
that = c\fi^, Qi = + Oi) and Pi = -^{a\ - ai). 

Thus the quantum Hamiltonian is 
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The coupling between the i*'' phononic mode and the 
/Lt*^ electronic state is given by 
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where oji is the frequency of the mode and is the en- 
ergy of the /J-*'' electronic state. Therefore, to calculate 
the electron-phonon coupling constants via the frozen 
phonon method one begins by calculating the electronic 
structure and the normal modes of the system; one then 
makes a small displacement of the system along each nor- 
mal coordinate, and thus calculates dX^/dQi. Therefore 
if there are N normal modes N such calculations must be 
performed. The frozen phonon method can be speeded 
up by the use of density functional perturbation theory 
(DFPT) 0. However, DFPT is an 0{N) technique, al- 
beit with a smaller coefficient than the frozen phonon 
calculation. Below we present a new method to calculate 
dXf_i/dQi in a single calculation. 
Janak's theorem 'lOl states that 
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where E is the total energy of the system and is the 
electronic occupancy of the state /z. It is important to 
realise that, within the density functional formalism 
is not required to be an integer It follows from Q 
and |(SJ| that 
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where {dE / dQi)x indicates that the derivative is taken 
after the charge is changed by x relative to the charge of 
the initially optimised geometry. We have used the fact 
that for a (meta)stable geometry {dE/dQi)o = 0. Once 
the electronic structure is solved in the equilibrium ge- 
ometry of the charge neutral system with a small change 
in the charge the forces on the molecule can be calcu- 
lated using the Hellman-Feynman theorem and, as the 
dynamical matrix is already known from the calculation 
of phonon spectrum, the electron-phonon coupling con- 
stants can readily be calculated. 
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FIG. 1: Comparison of the electron-phonon coupling con- 
stants for a variety of molecules calculated by the order-one 
method, g\, and the frozen phonon method g/ph (both in 
Ha/Bohr). It can be seen that the agreement between the 
frozen phonon and order-one methods is excellent. In the in- 
set we plot the same data on a log-log scale which displays 
the excellent agrement achieved even for small values of g. 
For each molecule we have calculated electron-phonon inter- 
actions associated with the HOMO and LUMO levels. 



Conceptually it is useful to realise that in the frozen 
phonon and DFPT approaches one considers the problem 
from the point of view of the electrons. That is one calcu- 
lates the electron-phonon coupling constants by making a 
small perturbation to nuclei and considering the resultant 
change on the electron. However, in the C(l) approach 
one asks the question what happens to the nuclei when 
one puts a small additional charge in the system. In most 
cases one is only interested in the coupling of the phonons 
to the electronic states closest to the Fermi level, how- 
ever, if one where interested in couplings to other states 
extensions of Janak's theorem to other electronic states 
allows for the determination of these parameters straight- 
forwardly within the 0(1) method. 

Thus the main result of this letter is that using equa- 
tion lO the gi can be calculated by solving a single nu- 
clear geometry. In the remainder of this work we present 
some benchmark calculations which show that equation 
© reproduces the values of the electron-phonon coupling 
constants calculated by the frozen phonon method for a 
number of small molecules. To do this we have imple- 
mented our scheme for calculating electron-phonon cou- 
pling constants in the Naval Research Laboratory Molec- 
ular Orbital Library (NRLMOL) El E E [11 H 
Il7| . Throughout we have used the Perdew-Burke- 
Ernzerhof (PBE) jlSj] exchange correlation functional. 
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FIG. 2; The ionisation potentials and electron affinities of 
the molecules considered in Fig. We plot the difference 
in energy, AE between the correct result (within the PBE 
approximation) and the results from approximating the ioni- 
sation potential (electron affinity) by estimating the Hubbard 
U and relaxation energy from the C(l) method. It can be seen 
that the approximation improves dramatically as the number 
of atoms in the molecule, TV, increases. The inset show the 
same data on a log-log scale, this suggests that the error in 
the approximation decreases as AE ~ iV"" with a ~ 1.5 — 2. 



In figure ^ we plot the calculated electron-phonon 
interactions for tetrathiafulvalene (TTF), tetra- 
cyanoquinodimcthane (TCNQ), bis(ethylenedithio)- 
tetrathiafulvalene (BEDT-TTF), Ceo and a number 
of diatomic and other smaU molecules calculated by 
the frozen phonon method against the electron-phonon 
coupling constants calculated by our method. It can be 
seen that the values calculated by either method are the 
same within numerical noise. 

Note that, as the normal coordinate is the eigenvector 
of the dynamical matrix, the sign of Qi is not well de- 
fined. Therefore, although dE is negative definite, the 
sign oi gi^ is also not well defined. However, the product 
QifidQi is well defined and describes the distortion caused 
by the addition of a charge. For example, for the stretch- 
ing mode in a simple diatomic molecule, the sign of the 
product gif^dQi indicates whether the bond is lengthened 
or shortened when a charge is added. 

Because the calculation of electron-phonon coupling 
constants is now computationally inexpensive, it is possi- 
ble to use the electron-phonon coupling to estimate other 
quantities. For example the adiabatic ionisation energy 
(or electron affinity) of a molecule, as opposed to the 
vertical ionisation energy (electron affinity) may be ex- 
tracted from the information already calculated. To do 
this one determines the Hubbard U (a second deriva- 
tive) from either the energy as a function of ^SUfj, and 
=F2(5n^ or from the appropriate eigenvalue as a function 



of ^SUf^ As the electron-phonon coupling constants are 
known the adiabatic ionisation energy (electron affinity) 
can be calculated directly from Markus-Hush theory. In 
figure|21we plot the difference between the adiabatic ioni- 
sation energies (electron affinities) calculated in this way 
and the correct result (within the PBE framework) for 
the molecules considered in Fig. as a function of the 
the number of the atoms in the molecules. While for 
small molecules this approximation is not accurate, the 
results for the larger molecules are excellent. Clearly it is 
for larger molecules that this approximation is needed as 
the computational power required to relax the geometry 
grows with the size of the molecule. 

It is important to stress that the calculation of the 
electron-phonon coupling constants is, in fact, intrinsi- 
cally an 0(1) problem as the dynamical matrix and the 
forces calculated on the geometries used to calculate the 
dynamical matrix together with the eigenvalue changes 
determined during the course of these calculations con- 
tain all the information required to calculate the gi^. 
However, significant computational complexities must be 
overcome to retrieve this information. Further, the calcu- 
lated electron-phonon coupling is rather sensitive to the 
size of the displacement used in a frozen-phonon calcu- 
lation. In general a simple criterion for the size of the 
displacement, e.g., x — y/2A/k where A is a small en- 
ergy, does not produce uniformly reliable results for all 
possible phonon frequencies in the frozen phonon calcu- 
lations of the type discussed above. However, we have 
found that the results for the 0(1) method are highly 
insensitive to the value of the 'small' charge, Sn^, used 
in the calculation. 

Finally we note that the calculation of isotope effects is 
an 0(0) problem in our method, i.e. only trivial matrix 
manipulation and no further solutions of the Schrodinger 
equation are required to calculate the electron-phonon 
coupling constants as the isotopic masses are varied. This 
contrasts to the frozen-phonon method where, as the dy- 
namical matrix and thus the normal modes change upon 
isotopic substitution, 0{N) additional calculations are 
required. 

A particularly interesting isotope effect is observed 
in the superconductor K-(BEDT-TTF)2Cu[N(CN)2]Br. 
When the eight hydrogen atoms in the BEDT-TTF 
molecule are replaced by deuterium the ground state is 
found to be a Mott insulator 19]. In the crystal the high- 
est occupied molecular orbital (HOMO) of the BEDT- 
TTF molecule is doped with holes from the anion layer, 
therefore one expects 3] that, in the small polaron limit, 
on deuteration, the bandwidth will change by a factor 
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where Wd (Wh) is the bandwidth of the deuterated (hy- 
drogenated) system and giho (gihH) is the coupling be- 
tween the i*^ mode and the HOMO in the deuterated 



3 



(hydrogenated) molecule. As both the deuterated and 
hydrogenated systems are very close to the Mott transi- 
tion [ijj we suggest that polarons may play a significant 
role in driving the Mott transition by dueteration. 

In conclusion we have demonstrated that by consider- 
ing the forces on the nuclei due to the addition or sub- 
traction of an arbitrarily small electronic charge one may 
calculate the elcctron-phonon coupling constants as an 
problem. This method is exact within the density 
functional formalism and was shown to be numerically 
accurate for a large number of small molecules. Note 
that, although we have only considered molecular sys- 
tems in our numerical work this is no intrinsic limitation 
of this method which prevents it being applied to infinite 
systems. 
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